3. Portfolio Theory
With Short Sale
library(quadprog)
R = 100*return[,2:13]
mean_p <- apply(R,2,mean)
cov_p <- cov(R)
sd_vect_p <- sqrt(diag(cov_p))
# min(mean_p)
# max(mean_p)
### With shortsale
M_p = length(mean_p)
Amat_p <- cbind(rep(1,M_p),mean_p)
mu_P = seq(0.07, 5.4, length = 300)
# Target portfolio means for the expect portfolio return
sd_P = mu_P # set up storage for std dev's of portfolio returns
weights_p = matrix(0, nrow = 300, ncol = M_p) # storage for return
for (i in 1:length(mu_P)) { # find the optimal portfolios
bvec_p <- c(1, mu_P[i])
result_p = solve.QP(Dmat = 2 * cov_p, dvec = rep(0, M_p), Amat = Amat_p,
bvec = bvec_p, meq = 2)
sd_P[i] = sqrt(result_p$value)
weights_p[i, ] = result_p$solution
}
plot(sd_P, mu_P, type = "l", xlim = c(0,15), ylim = c(0, 6), lty = 3, lwd = 2, main = "With Short Sale")
# plot efficient frontier (and inefficient portfolios below the min var portfolio)
mufree_p = mean(return$`Treasury Bill 3 month (rf)`)# input value of risk-free interest rate
points(0, mufree_p, cex = 4, pch = "*") # show risk-free asset
sharpe_p = (mu_P - mufree_p) / sd_P # compute Sharpes ratios
ind_p = (sharpe_p == max(sharpe_p)) # Find maximum Sharpes ratio
#weights_p[ind_p,] # print the weights of the tangency portfolio
lines(c(0, 15), mufree_p + c(0, 15) * (mu_P[ind_p] - mufree_p) / sd_P[ind_p], lwd = 4,
lty = 1, col = "blue") # show line of optimal portfolios
points(sd_P[ind_p], mu_P[ind_p], cex = 4, pch = "*") # tangency portfolio
ind2_p = (sd_P == min(sd_P)) # find the minimum variance portfolio
points(sd_P[ind2_p], mu_P[ind2_p], cex = 2, pch = "+") # min var portfolio
ind3_p = (mu_P > mu_P[ind2_p])
lines(sd_P[ind3_p], mu_P[ind3_p], type = "l", xlim = c(0, 25), ylim = c(0,30),
lwd = 3, col = 'red') # plot the efficient frontier
for(i in 1:12){
text(sd_vect_p[i], mean_p[i],colnames(return[,i+1]), cex=0.8)
}

### MVP
(mvp_meanreturn <- mu_P[ind2_p])
## [1] 1.816957
(mvp_sd <- sd_P[ind2_p])
## [1] 4.980327
weights_mvp <- weights_p[ind2_p,]
weights_mvp <- t(data.frame(weights_mvp))
colnames(weights_mvp) <- colnames(return[2:13])
knitr::kable(round(weights_mvp*100,2))
| weights_mvp |
50.17 |
-4.98 |
13.48 |
4.89 |
3.85 |
13.16 |
4.31 |
-9.16 |
21.82 |
8.07 |
-5.96 |
0.34 |
(mvp_meanreturn_ann <- mvp_meanreturn*12)
## [1] 21.80348
(mvp_sd_ann <- mvp_sd*sqrt(12))
## [1] 17.25236
### Efficient Portfolio Frontier
EPF_mean <- mu_P[ind3_p]
EPF_sd <- sd_P[ind3_p]
### Tangency Portfolio
(tan_meanreturn <- mu_P[ind_p])
## [1] 5.4
(tan_sd <- sd_P[ind_p])
## [1] 9.86121
(tan_var <- tan_sd^2)
## [1] 97.24347
(tan_sharpes <- (tan_meanreturn-mufree_p)/tan_sd)
## [1] 0.4753989
VaR&ES
### Tail dependence can be seen among the assets, therefore we can fit
### our portfolio with multivariate t-distribution.
## MVP VaR&ES
library(MASS)
alpha = 0.05
return_mvp <- rowSums(data.frame(
weights_mvp[1] * return[, 2],
weights_mvp[2] * return[, 3],
weights_mvp[3] * return[, 4],
weights_mvp[4] * return[, 5],
weights_mvp[5] * return[, 6],
weights_mvp[6] * return[, 7],
weights_mvp[7] * return[, 8],
weights_mvp[8] * return[, 9],
weights_mvp[9] * return[, 10],
weights_mvp[10] * return[, 11],
weights_mvp[11] * return[, 12],
weights_mvp[12] * return[, 13]))
fitt_mvp = fitdistr(return_mvp,"t")
param_mvp = as.numeric(fitt_mvp$estimate)
mean_mvpfit = param_mvp[1]
df_mvpfit = param_mvp[3]
sd_mvpfit = param_mvp[2] * sqrt((df_mvpfit) / (df_mvpfit - 2))
lambda_mvpfit = param_mvp[2]
qalpha_mvp = qt(alpha, df = df_mvpfit)
VaR_par_mvp = -100000 * (mean_mvpfit + lambda_mvpfit * qalpha_mvp)
es1_mvp = dt(qalpha_mvp, df = df_mvpfit) / (alpha)
es2_mvp=(df_mvpfit+qalpha_mvp^2)/(df_mvpfit-1)
es3_mvp=-mean_mvpfit+lambda_mvpfit*es1_mvp*es2_mvp
ES_par_mvp = 100000*es3_mvp
VaR_par_mvp
## [1] 6285.691
ES_par_mvp
## [1] 8957.941
## Asset VaR
S0 = 100000
qnalpha = qnorm(0.05)
### MSFT
q_msft = as.numeric(quantile(return$MSFT, alpha))
VAR_msft = -S0 * q_msft
#VAR_msft
### TSLA
fit_tsla <- fitdistr(return$TSLA, "normal")
param_tsla = as.numeric(fit_tsla$estimate)
mean_tsla = param_tsla[1]
sd_tsla = param_tsla[2]
VAR_tsla = -S0*(mean_tsla+qnalpha*sd_tsla)
#VAR_tsla
### AAPL
fit_aapl <- fitdistr(return$AAPL, "normal")
param_aapl = as.numeric(fit_aapl$estimate)
mean_aapl = param_aapl[1]
sd_aapl = param_aapl[2]
VAR_aapl = -S0*(mean_aapl+qnalpha*sd_aapl)
#VAR_aapl
### TWTR
fit_twtr <- fitdistr(return$TWTR, "normal")
param_twtr = as.numeric(fit_twtr$estimate)
mean_twtr = param_twtr[1]
sd_twtr = param_twtr[2]
VAR_twtr = -S0*(mean_twtr+qnalpha*sd_twtr)
#VAR_twtr
### AMZN
fit_amzn <- fitdistr(return$AMZN, "normal")
param_amzn = as.numeric(fit_amzn$estimate)
mean_amzn = param_amzn[1]
sd_amzn = param_amzn[2]
VAR_amzn = -S0*(mean_amzn+qnalpha*sd_amzn)
#VAR_amzn
### FB
q_fb = as.numeric(quantile(return$FB, alpha))
VAR_fb = -S0 * q_fb
#VAR_fb
### NFLX
q_nflx = as.numeric(quantile(return$NFLX, alpha))
VAR_nflx = -S0 * q_nflx
#VAR_nflx
### AAL
q_aal = as.numeric(quantile(return$AAL, alpha))
VAR_aal = -S0 * q_aal
#VAR_aal
### DAL
q_dal = as.numeric(quantile(return$DAL, alpha))
VAR_dal = -S0 * q_dal
#VAR_dal
### BAC
q_bac = as.numeric(quantile(return$BAC, alpha))
VAR_bac = -S0 * q_bac
#VAR_bac
### NVDA
q_nvda = as.numeric(quantile(return$NVDA, alpha))
VAR_nvda = -S0 * q_nvda
#VAR_nvda
### WBD
q_wbd = as.numeric(quantile(return$WBD, alpha))
VAR_wbd = -S0 * q_wbd
#VAR_wbd
VAR_asset <- c(VAR_msft, VAR_tsla, VAR_aapl, VAR_twtr, VAR_amzn, VAR_fb, VAR_nflx,
VAR_aal, VAR_dal, VAR_bac, VAR_nvda, VAR_wbd)
var_df <- cbind(names, VAR_asset)
var_df <- rbind(var_df,VaR_par_mvp)
var_df[13,1] <- c("MVP")
var_df[13,2] <- 6285.69101812641
colnames(var_df) <- c("Assets", "Var")
knitr::kable(var_df)
|
MSFT |
6918.17869693084 |
|
TSLA |
23502.7562034031 |
|
AAPL |
9898.52221650451 |
|
TWTR |
23279.5289061133 |
|
AMZN |
10763.0915501661 |
|
FB |
10561.0446590248 |
|
NFLX |
13072.710950181 |
|
AAL |
15706.0583956349 |
|
DAL |
12891.9240002285 |
|
BAC |
13041.1530377474 |
|
NVDA |
11974.7969532429 |
|
WBD |
14903.9099722124 |
| VaR_par_mvp |
MVP |
6285.69101812641 |
Assets’ Sharpe’s Ratios
names_sh <- data.frame(colnames(return[1,2:14]))
sharpes_list <- rep(NA, 13)
for(i in 2:14){
sharpes_list[i-1] = (mean(unlist(return[,i]))-mean(unlist(return[,15]))/100)/sd(unlist(return[,i]))
}
sharpes_list <- data.frame(sharpes_list)
shar_df <- cbind(names_sh,sharpes_list)
shar_df
## colnames.return.1..2.14.. sharpes_list
## 1 MSFT 0.30124000
## 2 TSLA 0.23644366
## 3 AAPL 0.26624907
## 4 TWTR -0.01416282
## 5 AMZN 0.23071594
## 6 FB 0.11372558
## 7 NFLX 0.15817121
## 8 AAL -0.05222523
## 9 DAL 0.01341556
## 10 BAC 0.08380025
## 11 NVDA 0.37645499
## 12 WBD -0.05347243
## 13 S&P500 0.08091924
Without shortshale
### Without shortsale
#R = 100*return[,2:13]
#mean_p <- apply(R,2,mean)
#cov_p <- cov(R)
#sd_vect_p <- sqrt(diag(cov_p))
### With shortsale
#M_p = length(mean_p)
Amat_p_noss <- cbind(rep(1,M_p),mean_p, diag(1,nrow=M_p))
mu_P_noss = seq(min(mean_p)+0.0001, max(mean_p)-0.0001, length = 300)
# Target portfolio means for the expect portfolio return
sd_P_noss = mu_P_noss # set up storage for std dev's of portfolio returns
weights_p_noss = matrix(0, nrow = 300, ncol = M_p) # storage for return
for (i in 1:length(mu_P_noss)) { # find the optimal portfolios
bvec_p_noss <- c(1, mu_P_noss[i], rep(0,M_p))
result_noss = solve.QP(Dmat = 2 * cov_p, dvec = rep(0, M_p), Amat = Amat_p_noss,
bvec = bvec_p_noss, meq = 2)
sd_P_noss[i] = sqrt(result_noss$value)
weights_p_noss[i, ] = result_noss$solution
}
plot(sd_P_noss, mu_P_noss, type = "l", lty = 3,
lwd = 2, xlim = c(0,15), ylim = c(0,7), main = "Without Short Sale")
# plot efficient frontier (and inefficient portfolios below the min var portfolio)
#mufree_p = mean(return$`Treasury Bill 3 month (rf)`) # input value of risk-free interest rate
points(0, mufree_p, cex = 4, pch = "*") # show risk-free asset
sharpe_p_noss = (mu_P_noss - mufree_p) / sd_P_noss # compute Sharpes ratios
ind_p_noss = (sharpe_p_noss == max(sharpe_p_noss)) # Find maximum Sharpes ratio
#weights_p[ind_p,] # print the weights of the tangency portfolio
lines(c(0, 15), mufree_p + c(0, 15) * (mu_P_noss[ind_p_noss] - mufree_p) / sd_P_noss[ind_p_noss], lwd = 4,
lty = 1, col = "blue") # show line of optimal portfolios
points(sd_P_noss[ind_p_noss], mu_P_noss[ind_p_noss], cex = 4, pch = "*") # tangency portfolio
ind2_p_noss = (sd_P_noss == min(sd_P_noss)) # find the minimum variance portfolio
points(sd_P_noss[ind2_p_noss], mu_P_noss[ind2_p_noss], cex = 2, pch = "+") # min var portfolio
ind3_p_noss = (mu_P_noss > mu_P_noss[ind2_p_noss])
lines(sd_P_noss[ind3_p_noss], mu_P_noss[ind3_p_noss], type = "l",
lwd = 3, col = 'red') # plot the efficient frontier
for(i in 1:12){
text(sd_vect_p[i], mean_p[i],colnames(return[,i+1]), cex=0.8)
}

### MVP
(mvp_meanreturn_noss <- mu_P_noss[ind2_p_noss])
## [1] 1.977063
(mvp_sd_noss <- sd_P_noss[ind2_p_noss])
## [1] 5.110456
weights_mvp_noss <- weights_p_noss[ind2_p_noss,]
weights_mvp_noss <- t(data.frame(weights_mvp_noss))
colnames(weights_mvp_noss) <- colnames(return[2:13])
weights_mvp_noss
## MSFT TSLA AAPL TWTR AMZN
## weights_mvp_noss 0.487966 -1.764314e-18 0.06592605 0.05609903 0.03247623
## FB NFLX AAL DAL BAC
## weights_mvp_noss 0.1346882 0.02256666 -2.931383e-17 0.146045 0.04453321
## NVDA WBD
## weights_mvp_noss 6.441445e-18 0.009699598
(mvp_meanreturn_ann_noss <- mvp_meanreturn_noss*12)
## [1] 23.72476
(mvp_sd_ann_noss <- mvp_sd_noss*sqrt(12))
## [1] 17.70314
sum(weights_mvp_noss)
## [1] 1
rownames(weights_mvp_noss) <- c("weight")
knitr::kable(round(weights_mvp_noss*100,4))
| weight |
48.7966 |
0 |
6.5926 |
5.6099 |
3.2476 |
13.4688 |
2.2567 |
0 |
14.6045 |
4.4533 |
0 |
0.97 |
### Efficient Portfolio Frontier
EPF_mean_noss <- mu_P_noss[ind3_p_noss]
EPF_sd_noss <- sd_P_noss[ind3_p_noss]
### Tangency Portfolio
(tan_meanreturn_noss <- mu_P_noss[ind_p_noss])
## [1] 3.999688
(tan_sd_noss <- sd_P_noss[ind_p_noss])
## [1] 7.96967
(tan_var_noss <- tan_sd_noss^2)
## [1] 63.51564
(tan_sharpes_noss <- (tan_meanreturn_noss-mufree_p)/tan_sd_noss)
## [1] 0.412526
alpha = 0.05
return_mvp_noss <- rowSums(data.frame(
weights_mvp_noss[1] * return[, 2],
weights_mvp_noss[2] * return[, 3],
weights_mvp_noss[3] * return[, 4],
weights_mvp_noss[4] * return[, 5],
weights_mvp_noss[5] * return[, 6],
weights_mvp_noss[6] * return[, 7],
weights_mvp_noss[7] * return[, 8],
weights_mvp_noss[8] * return[, 9],
weights_mvp_noss[9] * return[, 10],
weights_mvp_noss[10] * return[, 11],
weights_mvp_noss[11] * return[, 12],
weights_mvp_noss[12] * return[, 13]))
fitt_mvp_noss = fitdistr(return_mvp_noss,"t")
param_mvp_noss = as.numeric(fitt_mvp_noss$estimate)
mean_mvpfit_noss = param_mvp_noss[1]
df_mvpfit_noss = param_mvp_noss[3]
sd_mvpfit_noss = param_mvp_noss[2] * sqrt((df_mvpfit_noss) / (df_mvpfit_noss - 2))
lambda_mvpfit_noss = param_mvp_noss[2]
qalpha_mvp_noss = qt(alpha, df = df_mvpfit_noss)
VaR_par_mvp_noss = -100000 * (mean_mvpfit_noss + lambda_mvpfit_noss * qalpha_mvp_noss)
es1_mvp_noss = dt(qalpha_mvp_noss, df = df_mvpfit_noss) / (alpha)
es2_mvp_noss=(df_mvpfit_noss+qalpha_mvp_noss^2)/(df_mvpfit_noss-1)
es3_mvp_noss=-mean_mvpfit_noss+lambda_mvpfit_noss*es1_mvp_noss*es2_mvp_noss
ES_par_mvp_noss = 100000*es3_mvp_noss
VaR_par_mvp_noss
## [1] 6291.464
ES_par_mvp_noss
## [1] 9034.092
##4. Asset Allocation
Efficient Portfolio
#efficient portfolio
new_dat = return[,2:13]
return_year = 12*colMeans(new_dat)
cov_mat = cov(new_dat)
eff_port = efficient.portfolio(er = return_year,cov_mat,target.return = 0.06,shorts = FALSE)
eff_port$er/12
## [1] 0.004999994
#check return
ret=0
for (i in 1:12)
{
val = as.numeric(return_year[i])*as.numeric(eff_port$weights[i])
ret = ret+val
}
ret
## [1] 0.05999993
#monthly risk
eff_risk = eff_port$sd/sqrt(12)
eff_risk
## [1] 0.02134934
# vaR = -S0*(u + qnorm(0.05)*sd)
S0 = 100000
q = qnorm(0.05)
eff_var = -S0*(eff_port$er/12 + q*eff_risk)
eff_var
## [1] 3011.654
#excel
cut_off = qnorm(0.05,mean = (1+eff_port$er/12)*S0,sd = S0*eff_risk)
var2 = S0-cut_off
var2
## [1] 3011.654
var2/S0
## [1] 0.03011654
#Expected Shortfall
S0*((-eff_port$er/12) + (dnorm(qnorm(0.05))/0.05)* eff_risk)
## [1] 3903.755
Tangency portfolio
risk_free = colMeans(return[,15]/100)
tan_port = tangency.portfolio(er=return_year, cov.mat=cov_mat, risk.free=risk_free,shorts=FALSE)
tan_port
## Call:
## tangency.portfolio(er = return_year, cov.mat = cov_mat, risk.free = risk_free,
## shorts = FALSE)
##
## Portfolio expected return: 0.4214667
## Portfolio standard deviation: 0.06885209
## Portfolio weights:
## MSFT TSLA AAPL TWTR AMZN FB NFLX AAL DAL BAC NVDA
## 0.4830 0.0686 0.0903 0.0000 0.0452 0.0000 0.0000 0.0000 0.0000 0.0000 0.3129
## WBD
## 0.0000
# w1 for tangency portfolio, w1 for risk free
# w1+w2 = 1
# 0.4799845 * w1 + 0.08543894 * w2 = 0.06
# 1.3041868% for the tangency portfolio and -0.3041868% for the risk free
A = matrix(c(1, tan_port$er, 1, as.numeric(risk_free)), ncol=2)
B = matrix(c(1, 0.06),ncol = 1)
solve(A)%*%B
## [,1]
## [1,] 0.1276228
## [2,] 0.8723772
#vaR
S0 = 100000
q = qnorm(0.05)
eff_var = -S0*(0.05/12 + q*0.06885209/sqrt(12))
eff_var
## [1] 2852.626
eff_var/S0
## [1] 0.02852626
#ES
S0*((-0.05/12) + (dnorm(qnorm(0.05))/0.05)* 0.06885209/sqrt(12))
## [1] 3683.158
5. PCA
correlation matrix of the returns on your assets.
AAL and DAL are the assets that most highly correlated with
correlation = 0.78. If we only consider about the 12 assets, then NFLX
and DAL has the lowest correlation = 0.0089. The NFLX also has a
correlation with Treasury Bill 3 month (risk free) = 0.006.
Diversification will reduce risk for the assets that are negatively
correlated. For example, if the correlation between 2 assets is -1,
which is the perfect value for diversification, if one asset goes up 5%,
the other asset will goes down 5%. So combining these 2 assets would
minimize our risk to 0. In our 12 assets, we did not see any negative
correlation. As a result, diversification will not reduce risk.
(note: here are a few assets that are slightly correlated such as
NFLX and DAL, NFLX and AAL, TWTR and AAL, NVDA and WBD, MSFT and
TWTR.)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:psych':
##
## logit
## The following object is masked from 'package:fBasics':
##
## densityPlot
df <- return[,2:13]
# correlation matrix
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
usr <- par("usr")
on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
Cor <- abs(cor(x, y))
txt <- paste0(prefix, format(c(Cor, 0.123456789), digits = digits)[1])
if(missing(cex.cor)) {
cex.cor <- 0.4 / strwidth(txt)
}
text(0.5, 0.5, txt,
cex = 1 + cex.cor * Cor)
}
pairs(df,
upper.panel = panel.cor, # Correlation panel
lower.panel = panel.smooth) # Smoothed regression lines
Parallel coordinate plot
parcoords(df, rownames = F, reorderable = T , queue = T,
withD3 = TRUE, alpha=0.5)
PCA Analysis
Another way to dig into the correlation between the return of the
assets, we looked at PCA and its plots. PCA helps to reduce the
dimension of the dataset and compute principal components that contain
most of the information in the data. From the PCA Analysis, we found
that first eight components captures almost 90% of the total variance,
while nine components could explain 93% of the variance. The Scree Plot
shows that the unexplained part drop quickly by adding PC1 to PC6, and
gradually slower after that. Therefore, using around 6 to 9 PCs would
explain the majority of the variance in our data.
To further developed the correlations, we also draw a biplot.
According to the principle of the biplots, the variables are positively
correlated if their vectors form an angle smaller than 90 degrees. From
our plot, all the assets are having a acute angle, so we conclude that
they all having a positively relationship.
df = return[,2:13]
model_pca = prcomp(df, center = TRUE, scale. = TRUE)
su = summary(model_pca)
su
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.0761 1.3586 1.06307 0.93727 0.8879 0.80884 0.78736
## Proportion of Variance 0.3592 0.1538 0.09418 0.07321 0.0657 0.05452 0.05166
## Cumulative Proportion 0.3592 0.5130 0.60717 0.68038 0.7461 0.80060 0.85227
## PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.70899 0.67221 0.56977 0.54539 0.44294
## Proportion of Variance 0.04189 0.03766 0.02705 0.02479 0.01635
## Cumulative Proportion 0.89415 0.93181 0.95886 0.98365 1.00000
screeplot(model_pca, type = "line", main = "Scree plot")

autoplot(model_pca, loadings = TRUE, loadings.colour = 'blue',
loadings.label = TRUE)

#model_pca$rotation[,1:2]
factor analysis
X = scale(df,center=T,scale=T)
n=fa.parallel(X)

## Parallel analysis suggests that the number of factors = 2 and the number of components = 2
df.fa = fa(X,nfactors = 4,fm="mle",rotate = "varimax")
fa.diagram(df.fa)

faa1 = factanal(df,factor=3,rotation = "varimax")
faa1
##
## Call:
## factanal(x = df, factors = 3, rotation = "varimax")
##
## Uniquenesses:
## MSFT TSLA AAPL TWTR AMZN FB NFLX AAL DAL BAC NVDA WBD
## 0.459 0.744 0.534 0.853 0.411 0.617 0.626 0.230 0.193 0.502 0.553 0.338
##
## Loadings:
## Factor1 Factor2 Factor3
## MSFT 0.673 0.208 0.210
## TSLA 0.462 0.190
## AAPL 0.639 0.201 0.130
## TWTR 0.148 0.348
## AMZN 0.752 0.150
## FB 0.568 0.237
## NFLX 0.598 0.122
## AAL 0.148 0.855 0.130
## DAL 0.880 0.171
## BAC 0.298 0.488 0.414
## NVDA 0.634 0.201
## WBD 0.101 0.220 0.776
##
## Factor1 Factor2 Factor3
## SS loadings 2.870 1.934 1.135
## Proportion Var 0.239 0.161 0.095
## Cumulative Var 0.239 0.400 0.495
##
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 40.08 on 33 degrees of freedom.
## The p-value is 0.185